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Abstract —The empirical mode decomposition (EMD) has 
achieved its reputation hy providing a mniti-scale time-frequency 
representation of nonlinear and/or nonstationary signals. To 
extend this method to vector-valued signals (VvS) in multi¬ 
dimensional (multi-D) space, a multivariate EMD (MEMD) has 
been designed recently, which employs an ensemble projection to 
extract local extremum locations (LELs) of the given VvS with 
respect to different projection directions. This idea successfully 
overcomes the problems of locally defining extrema of VvS. 
Different from the MEMD, where vector-valued envelopes (VvEs) 
are Interpolated based on LELs extracted from the 1-D projected 
signal, the vector-valued EMD (VEMD) proposed in this paper 
employs a novel back projection method to interpolate the VvEs 
from 1-D envelopes in the projected space. Considering typical 
4-D coordinates (3-D location and time), we show by numerical 
simulations that the VEMD outperforms state-of-art methods. 

Index Terms —Empirical mode decomposition, intrinsic mode 
function, vector-valued signal, back projection, optimization. 

I. Introduction 

The empirical mode decomposition (EMD) was firstly de¬ 
signed for nonlinear and/or nonstationary signal analysis 111. 
Combined with the Hilbert transform, the Hilbert-Huang- 
Transform (HHT) 111 provides a finer time-frequency spectrum 
of a given signal compared with other well-known methods 
such as the Eourier transform, the wavelet transform, or 
the Wigner-Ville transform. Moreover, the EMD does not 
require any pre-defined basis. It decomposes a given signal 
/(/) into a finite number of intrinsic mode functions (IMEs), 
fj{t) := aj{t) cos{9j{t)), j = 1,..., J, and a monotonic trend 
rj+iit), i.e., 

J 

f{t):=J2fdt) + rj+i{t), tG[0,T]. (1) 

/=! 

Here, each fj{t) might be considered as an amplitude- 
modulated and frequency-modulated signal, or as a generalized 
Eourier component 111. The properties of the signal model 
are extensively studied in ID, 11 and the references therein. 
Based on the instantaneous frequency of each IME is 
well defined by ojjit) := Q. In addition, the corre¬ 

sponding Hilbert amplitude spectrum can be naturally derived 
as {aj{t) on the curves t S [0,T], j = 1,..., J} 

0. Up to now, the EMD and its many variations have been 
successfully employed in many disciplines, such as signal 
processing El, hydrology 16), and geophysics [|7]. 
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In general, the EMD aims at sequentially extracting each 
IME fj{t) through a filtering operation called sifting process. 
Considering the signal at the /cth level, Xk{t) := f{t) - 
the sifting operator S can recursively be defined 
by 5"[a:fe](f) := where n is 

the iteration number, := Xk(t), and Al[-] represents 

the local trend approximation operator. In the classic EMD 0, 
M. [a:] (f) is defined as the mean curve of the upper and lower 
envelopes which are defined by cubic spline interpolation of 
the local maxima and local minima of x{t), respectively. The 
sifting process stops at n = N until some ad hoc criterion 
is met, e.g. fk{t) := mimics in some sense a 

generalized Eourier component. 

The EMD is a completely data-driven decomposition 
method which heavily depends on the definition of local 
extrema. Since their definition is unclear in higher dimensions, 
or the extrema are nonunique, it is, therefore, difficult to 
extend the method to these cases. The original signal or 
univariate time series decomposed in Q consisted of data of 
the form / : [0, T] —K, mapping one-dimensional data in 
one-dimensional space. In higher dimensions, we distinguish 
between multivariate and vector-valued data. The multivariate 
case consists of ‘cube’ data, e.g., an image in 2-D, or a 
volume in 3-D. In this case, the function is of the form 
/ : [OjT]"^ -A K with d = 2,3 being the space dimension. 
Here, a local extremum may be defined as the strict extremum 
in a pre-defined neighborhood 0. However, in the vector¬ 
valued case, the data are multi-D ‘curve’ data which are of 
the form / : [0, T] -A In this case, the definition of local 
extrema is much more complicated since even the notion of 
‘neighborhood’ is unclear. 

In this paper, we concentrate on the latter, more difficult, 
case. We define the vector-valued signal as follows. 

Definition 1.1: Given a finite number of 1-D signals 
S C(M.) f] Lao{M.),i = l,...,d, t G [0,T], the cor¬ 
responding vector-valued signal (VvS) is defined as F{t) := 
i.e., F:[0,r]^K‘^. 

To develope an EMD-like decomposition for VvS, a 
straightforward idea is to decompose a complex signal, 
F{t) := by applying the classic EMD 

to the real and imaginary parts separately. However, this 
usually leads to a different number of IMEs for each of 
the two components which is an undesired effect a. In 
fact, if we reconsider the signal decomposition model in 3- 


D (2-D location and time), and then define the correspond¬ 
ing jth complex IMF as Fj{t) := = 

cos{0^^\t)), cos{9^^\t))'^, we will observe an 

interesting fact that each Fj (t) is a rotation invariant compo¬ 
nent. Its 3-D envelope {a^j^ (^)j 0 should be some tube¬ 

shaped surface tightly enclosing the 3-D curve Fj{t). Then 
its 3-D local trend/mean J\4[Fj]{t) can be considered as the 
barycenter curve of the 3-D envelope surface. This observation 
implies a possible way to approximate A4[Fj]{t) from the 
view of statistics. This means that we should interpolate 
the envelopes based on the local extrema of Fj{t) along a 
selected projection direction, and then average all interpolated 
envelopes with respect to all possible projection directions, 
following a idea from ifTOll . Here, we define within a general 
concept the multi-D local extremum as follows. 

Definition 1.2: Given a VvS F{t) and a selected unit 
projection direction p, the local extremum of Fit) along 
p is defined as the hyperpoint {F(tk),tk), where tk is the 
corresponding local extremum location (LET) of the projected 
1-D signal Vp[F]{t). 

If the projection number approaches infinity, M.[Fj\{t) 
should be the expectation of the local mean approximation 
based on the projection. Such operations belong to a typical 
ensemble approach. This implies that the uniform direction 
sampling scheme should be an optimal choice ini, ina. 
Moreover, decomposition tests with white noise show that such 
EMD-like decomposition of VvS perfectly inherit the dyadic 
filter bank property of the classic EMD lfT3l , lfT4ll . 

To further study the above mentioned local mean ap¬ 
proximation, we consider the vector-valued envelope (VvE) 
interpolation in a different but more general way than in 
Eol, m, M- Here, the methods directly interpolate the 
VvE in the multi-D space and assume that the following 
property is satisfied; if we project the interpolated VvE along 
the corresponding projection direction, the projected curve 
should be nothing but the envelope interpolated based on the 
local extrema in the projected 1-D space. Unfortunately, this 
assumption only meets the requirement of the direct inter¬ 
polation method (without constraints), e.g., the cubic spline 
interpolation (CSI), but might be defective for the others, e.g. 
optimization based method a, when additional constraints 
have to be maintained. Thus, a more natural way to obtain the 
VvE is to firstly interpolate the 1-D envelope in the projected 
space with the desired method, and then back-project it into 
the original multi-D space. 

In this paper, we explain how to realize the vector-valued 
EMD (VEMD) with an optimization based back projection. 
To simplify the discussion, a typical signal decomposition 
model in 4-D space (3-D location and time) will be studied. In 
the following. Section II recalls the recently developed multi¬ 
variate EMD (MEMD) with approximately uniform direction 
sampling. Section III introduces the proposed VEMD, Section 
IV presents the numerical studies, and Section V concludes 
the paper. 


H. Multivariate empirical mode decomposition 

To keep the terminology from previous papers, in this 
section, the MEMD concerns the VvS decomposition problem. 
Eor a concrete 4-D decomposition problem, the signal model 
0 can be generalized as 


Fit) :=J2F,it) + Rj+,it), 
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where F{t) := (f^^^t),/^“^^t), f^^\t))"’" is the given 4-D 
curve, Fj{t) := is the jth decom¬ 
posed IME, and Rj+i{t) := is 

the monotonic trend. 


A. Approximately uniform sampling on a unit sphere 

To obtain equidistributed direction on a unit sphere, in this 
paper, we employ a low-discrepancy sampling scheme based 
on transformed Hammersley points, which can produce more 
uniform samples on the hypersphere than other methods, like 
the polar coordinate lattices or rotation method ifTSll . 

Simply speaking, for any assumed direction number M and 
any prime base b,b > 2,b G N, a selected direction index 
m,m G {0,1, • • • has a unique digit expansion 

r-m 

m = '^CjnjF, (3) 

3=0 

where each Cm, is an integer in [0,5— 1]. The Van der Corput 
sequence can be defined as 


{zbim)\zbim) :='^Cmjb ^ \ m = 0,1, ■ • ■ , M - 1}, (4) 

j=o 

and the corresponding Hammersley points set is m 

,m = 0,1,-- - (5) 

Since G [0,1), and Zb{m) G [0,1), we can transform 
the equidistributed square [0,1) x [0,1) into equidistributed 
cylinder [0,27r) x [—1,1) by setting cj) := 2 tt^ and z := 
2zb{m) — 1. Einally the uniformly sampled projection direc¬ 
tions are pm '■= (vT2T^cos((/)), vT^T^) sin(())), m = 
0,1, • • • , M. Eig. shows two sampling examples for M = 
512, and 5 = 2,5 respectively. 

It should be noted that a more satisfactory sampling scheme 
in multi-D space can be generated by introducing different 
prime bases which connects to the Hammersley points based 
on Halton sequences na. 






Fig. 1. Two sampling examples based on Hammersley points, (a) 
Pm e (b) Pm e 


Algorithm 1 : MEMD in 4-D space 

1: set R\(t) = F{t), j = 1, and the values of M and b\ generate the 
projection direction set 

2: extract j'th IMF by sifting process {Fj(t) := 

(A) set intermediate curve F[jfi{t) = Fij{t) and n = 0; 

(B) approximate the local mean := 

(a) for every direction Pm £ = 0,1, ■ • ■ ,M — 1, project 

into 1-D space, ~ Fpm[Hj,n]{t)', 

(b) detect all LELs of (t) (maximum location and 

minimum location , fc“ £ N) ; 

(c) generate upper [/J™ (t) and lower (t) VvEs by 

interpolating data’ pairs 

k- )’ n k- ) CSI; ’ ’ 

(d) set TUt) ■■= ^ tmZo Uf-it) + Vf- it)- 

(C) update JTj,„+i(f) = Hj^nit) — Tj^n{t) and n = n + 1; 

(D) calculate stopping criterion SDj,„ as the one in (HI; 

(E) repeat steps (B) to (E) until SDj,„ < SDihr, or the max iteration 
number N is met; define j\h IMF Fj(f) := 

3: update = Rjit) — Fj{t), and j = j -\- 1; 

4: repeat steps 2 to 3 until number of extrema in Rj (t) is less than 
2 or an expected IMF index is met, i.e., j = J. 


B. Multivariate empirical mode decomposition 

The recent MEMD method has been shown its strength for 
VvS decomposition, especially for mode alignment problem 
and noise-assisted applications. The white noise decomposi¬ 
tion test illustrates its remarkable dyadic filter bank property 
comparing to classic EMD and ensemble EMD m. The 
whole algorithm can be summarized in Algorithm [T] 

In MEMD, the local mean approximation Jvi [•] (f) could be 
considered as an ensemble approach, which means the idea 
mean curve T’(f) may be exactly approximated as M —>■ oo. 
In other words, we could re-define the local mean by T'(f) 
M.[H]it) := E(i([/P™(f) + yP’^it)), where E denotes the 
expectation over the direction pm- In real application, the 
projection number M needs not to be a large number com¬ 
promising both computational complexity and approximation 
performance. 

III. Vector-valued empirical mode decomposition 

In this section, a different EMD extension, namely the 
VEMD, will be explained, in which the VvEs of the given 
VvS are generated by back-projecting the 1-D envelopes inter¬ 
polated in the projected space to the original 4-D space. Since 


the naive back-projection might result in infinite solutions, a 
novel optimization scheme will be designed to guarantee an 
unique VvE which should be also as smooth as possible. 

A. Optimization based back projection 

Keeping the notation system in algorithm [T] but ignoring 
the trivial subindexes, in this subsection, H{t),pm denote the 
considering VvS and the projection direction, Up^{t)^Vp^{t) 
denote the upper and lower VvEs w.r.t pm in the original 
4-D space, and Up^ (t), Vp^ (t) denote the upper and lower 
envelopes w.r.t pm in the 1-D projected space. 

Assuming H (t) and pm are given, and the Vp^ [H] (t) is the 
projected 1-D signal, the local maxima/minima of Vp^[H]it) 
can be easily detected at the corresponding time locations 
and tk-, ,k~ £ N. With some interpolation method, e.g. 
the CSI, the 1-D upper Up^{t) and lower Vp^it) envelopes 
can be interpolated based on the pairs iVp^[H]itk+),tk+) and 
i'Ppm tk-) separately. Now, we aim to project Up^ (f) 

and Vp^ (t) back to the original 4-D space in order to obtain 
the VvEs Up^{t) ■= and V(f) := 

w.r.t the direction pm- 

Mathematically speaking, such back projection problem 
is equivalent to the solution problem of the linear system 

'PpmPpmW) = '‘J-Pmit) (O^ 'Ppmi^PmW) = WWch 

has infinite solutions because the number of unknowns is 
larger than the number of equations. However, the infinite 
solutions can be constrained to unique one if next two 
properties can be maintained simultaneously; a) the back- 
projected VvEs Up^{t) (or Vp^{t)) should pass through the 
LELs iHitk+),tk+) and w.r.t direction pm', 

b) the VvEs Up^ (t) (or Vp^ (t)) should be as smooth as 
possible. To investigate the smoothness of an unknown 4-D 
VvS X{t) := we may employ an nth 

order Sobolev norm funtional, i.e., the L 2 norm of the nth 
(weak) derivative of the function X 

5M[X1 := |pW[xW]||i^ + ||PW[a;[2)]||2^ + ||PW[*[3]]||2^, (6) 

where is the nth order derivative operator as its matrix 
form for discrete problem is well-known . 

Now, the back-projected VvEs can be uniquely obtained by 
solving the following optimization problems 
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In fact, (PI) and (P2) are both quadratic optimization problems 
with equality constraints. They can be written as quadratic 
optimization problems without constraints by solving the 
linear constraint system and then implementing the variable 
reduction. 













B. Vector-valued empirical mode decomposition 


Comparing to the MEMD method in section II-B the 


VEMD employs the back projection to obtain the VvEs from 
the envelopes interpolated in the 1-D projected space w.r.t a 
selected projection direction. 

Most of the computation steps in VEMD algorithm are as 
the same as the ones in algorithm [T] except the step 2 :(B);(c) 
which should be replaced by 


Algorithm 2 ; VEMD in 4-D space 

2: (B):(c) generate upper and lower u?™ (i) envelope 

in the 1-D projected space by interpolating data pairs 

and 

based on CSI; then solve the corresponding back projection opti¬ 
mization problem (PI) and (P2) to obtain the 4-D upper (t) 
and lower it) envelopes; 


It should be noted that in algorithm the envelope inter¬ 
polation in the 1-D projected space can be implemented by 
any reasonable interpolation method in order to meet different 
mathematical requirements, e.g. the method in ||4|. 

IV. Numerical studies 

In this section, all the simulations are implemented in 
MATLAB on a laptop equipped with an i7-4700 quad-core 
CPU, 8 GB memory and under Windows. All optimization 
problems are solved by the standard CVX toolbox from M- 
To distinguish the associated variables in MEMD and VEMD, 
we employ the subscripts M and V. 

In the VEMD, the order of the derivative operator in (j^, 
the projection direction number and the prim base in Q are 
free parameters, each of which may effect the behavior of the 
method. To determine each one, we study the decomposition 
problem of the following VvS 

F{t) ■.= X(t) + Y{t), te [0,1], 

/ \ / (1-I-cos(27rt)) sin(207rt) \ 

X{t) := I I = I (2-I-cos(47rt)) sin(207rt) | , 

\ ) \ (3-|-cos(67rt))sin(207rt) j (7) 

/ \ / sin(47rt) \ 

Y{t)-.= ^ 2sin(47rt) 

V / \ 3sin(47rt) j 

Comparing to the signal model in (j^, the components X{t) 
and Y{t) can be considered as the hrst IMF Fi{t) and 
the residual R 2 {t), respectively. The reason that we set a 
common frequency of each a;!®! (t) and (f) is to facilitate the 
discussion. Since the MEMD, as the same as the VEMD, has 
the mode alignment property (see cax different frequencies 
involved in the components a;W(f) or undoubtedly 

lead to many more decomposed IMFs, each of which should 
contain a particular frequency. 

A. Derivative order 

Selecting a unit projection direction p = [^, |, ^], 

we can easily obtain the projected VvS 'Pp[F]{t) to¬ 
gether with its LELs (Figj^ (d)). Based on these 




Time t/s (e) 


Eig. 2. The envelope interpolation in MEMD. (a)-(c) The VvS F{t) 
and its interpolated VvEs Uuit) and 141(1:) w.r.t an unit projection 
direction p; (d) The projected signal 'Pp[F]{t) and its interpolated 
envelopes Up{t) and Vp{t); (e) The absolute difference between 
the projected VvEs and the envelopes in the projected space, e.g. 
Dum(1) := TpIUmKI) - Up{t)\- 


locations, by using a CSI, we can interpolate the en¬ 
velopes Up{t),Vp{t) in the projected space (Figj^ (d)). 


and the VvEs (7 m(1) := 




and 




M V'-U “M 

in the original 4-D 


space (Fig|^ (a)-(c)). Fig|^ (e) illustrates that the projected 
VvEs {V\UM]{t),F\VM]{t)) are nothing but the envelopes 
{up{t), Vp{t)) in the projected space. In other words, Figj^ 
graphically explains why the MEMD interpolates the VvEs 
directly in the original space but not in the projected space. 

Now, based on the interpolated envelopes together with the 
LELs in the projected space, we can obtain the 4-D VvEs 
using the proposed back projection method, i.e. solving the 


optimization problems (PI) and (P2) described in Sec III-A 


Considering that, normally, the first order derivative won’t 
be an optimum choice for smoothness measurement IH, we 
set the 2nd or 3rd order derivative in •ID alternatively for 
simulation. Figj^ presents the corresponding solved VvEs 
based on the data shown in Fig|^ (d). 

In Fig|^ sub-figures (a)(c)(e) illustrate back-projected VvEs 
U\/^(t)^Vy^{t) and where the subscripts V 2 

and V 3 denote the 2nd and 3rd order derivative selected 
in VEMD each. To evaluate the interpolation performance, 
we take the VvEs interpolated in MEMD as a benchmark. 
Sub-figures (b)(d)(f) present the absolute differences between 
each VvE in VEMD and the corresponding VvE in MEMD, 
e.g. Dt/MV 2 (i) := (D: 


''MV2 






lUuit) — ( 7 v 2 (f)|- These figures imply that the VEMD with 
2nd order derivative would be much close to the MEMD 
comparing to the VEMD with 3rd order derivative. This 
interesting phenomenon can be effortlessly understood if we 
can recall that the CSI requires the interpolated curve to be 
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Fig. 3. The envelope interpolation in VEMD. Left column (a)(c)(e): 
the VvS and its VvEs obtained by solving back projection problems 
(PI) and (P2) w.r.t an unit projection direction p. Subscripts V 2 and 
V 3 of the solution denote the 2nd and 3rd order derivative selected in 
for VvE interpolation; Right column (b)(d)(f): absolute difference 
between interpolated VvEs in VEMD and the ones in MEMD shown 
in Fig[2] e.g. DITmvsF) “ lUuit) - f/vaF)!; Last row (g)(h): The 
absolute difference between the projected VvEs and the envelope in 
the projected space, e.g. D«v 2 (f) := |Tp[( 7 v 2 ](f) ~ tip(f)|. 


at most 2nd order continuously differentiable. Therefore, the 
VEMD with 2nd order derivative can produce similar results 
to the ones from MEMD, while the VEMD with 3rd order 
derivative can provide smoother interpolated curve. Einally, 
sub-figures (g)(h) imply both VvEs interpolated in VEMD 
satisfy the first equality constraint in (PI) and (P2) with high 
accuracy. 


B. Projection number and prime base 


In section II-B we have shown that, given a fixed prime 
base for projection sampling, the local mean approximation 
operator M[-] might be considered as an ensemble approach. 
In practice, we have to determine a finite projection number 
such that the approximation would be good enough or at 
least approximately convergent. Let’s consider again the VvS 
F(t) in (j^. If we apply the operator on F[t) with a fixed 
projection number M, M.^\F]{t) shall be the approximated 
local mean, e.g. M^[F]{t), the result of step (B) in the 
MEMD (algorithm or the M.y^[F]{t) or in 

the VEMD (algorithm |^. To measure the approximation 
efficiency, we define a percent root mean squared difference 
(PRD) as an error function of the projection number M, 


PRD(M) := 


\\ymF 


X 100%. (8) 


SS5 

Q 

^4 


^PRDm(M) oPRDv2(M) •PRDv3(M) 




|*PRDm(6) oPRDv2(&) -PRDvsW 






5 250 500 750 . ,1000 

Direction Number M 


16 32 48 ,,,64 

Prime Number b 1'^.' 


Fig. 4. Local mean approximation performance of operators 
Af ^ [■] and M^l [•] in MEMD and VEMD. (a) PRD values of the 
three operators for b — 2, and M G [5,1000]; (b) PRD values of the 
three operators for M = 512, and b £ [2,80]. 


Pig|^ (a) describes the approximation performance of the 
operators Ad^[-] and Al^[-] for 5 = 2,M S 

[5,1000]. As can be seen, all PRD sequences decay dramati¬ 
cally in the beginning, and each one might be approximately 
convergent when M > 256. Since the sampled projection 
directions are not ideally uniform distributed on the sphere 
ca, the sampling error may result in slight oscillations on 
the PRD sequence. Therefore, we set M = 512 for rest 
simulations. 

Pig|^ (b) presents the performance of all three approxima¬ 
tion operators with different prime bases b G [2,80]. It shows 
that different prime numbers do not effect the PRD values 
significantly. This is because the prime base only effects on 
the projection locations on the sphere (see Fig{^. When the 
projection number is large enough, approximation operator 
with different prime bases should provide consistent perfor¬ 
mance. In other words, the prime base can not be set as some 
number which is close to the projection number. Otherwise, 
from © and (|^, we will find that Zi,{m) and m may have 
strong correlation, e.g. partial linear dependence, that results 
in a non-uniform sampling on the sphere. Therefore, we set 
6 = 2 for rest simulations. 

Fig|^ illustrates another exciting phenomenon; the approx¬ 
imation operator works better than other two for 

arbitrarily selected projection number and prime base. This 
can be easily understood with following two considerations: 
1 ) Al^[-] requires higher order derivatives for smoothness 
evaluations than other two; 2) many common functions are 
continuously differentiable with high order. On the other hand, 
the performance of [•] and the one of [•] coincide 
each other perfectly, which further implies the effectiveness 
of the proposed method. 

C. Decomposition performance 

With the fixed parameters, the proposed VEMD can be 
applied to decompose the given VvS F(t) in 0. The corre¬ 
sponding decomposition results are shown in Eig|^ As can 
be seen, all decomposed components by using MEMD or 
VEMD with different derivative orders are very similar to 
the ideal ones. To numerically distinguish the decomposition 
performance, the PRDs of the decomposition results w.r.t each 
method are listed in Table-|I] The numerical results support 
our observation again: the VEMD with 3rd derivative order 
is better than the one with 2nd derivative order and the 
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Fig. 5. Decomposition results of the given VvS F(t) in ([^ by using 
MEMD and VEMD methods. Left column (a)(c)(e): the original com¬ 
ponent X{t), the first decomposed Imfs using MEMD (Fi^uit)), and 
VEMD with 2nd derivative (Fiy^{ty) and 3rd derivative (Fi,v 3 (f)), 
and the absolute difference between X{t) and each Fi^.{t), e.g. 
DE’i,M(f) := \X{t) — fi,M(i)|; Right column (b)(d)(f): the original 
component Y{t), the residual i? 2 ,M(i), f? 2 ,V 2 (f) and R 2 y^{t), and 
the corresponding absolute difference. 


TABLE I 

PRD OF THE DECOMPOSITION RESULTS USING MEMD AND VEMD 


PRD(%) 

M 

V2 

V3 

FiAt) 

2.95 

2.95 

2.65 

Ra.-W 

3.10 

3.10 

2.79 


MEMD; and the latter two methods have similar behavior from 
mathematic point of view. 

V. Conclusion 

As a representative data-driven method, EMD can provide 
finer time-frequency analysis of any nonlinear and/or non¬ 
stationary signal sampled in (non-)uniform grids. This paper 
introduced a novel method to extend the classical EMD for 
vector-valued signal decomposition. Different from the exist¬ 
ing MEMD, our proposed VEMD obtains the local mean curve 
by projecting the envelopes obtained in 1-D projected space 
back into the original multi-D space. Since the VEMD does 
not require any particular envelope interpolation method, it is 
more general and flexible compared to the MEMD for many 
potential applications. In addition, with the method introduced 
in HU, we can generate a time-frequency representation of 
any VvS which can then provide meaningful information for 
further vector-valued data analysis. 
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